library(readxl)
library(ggvenn)
library(clusterProfiler)
source('00_util_scripts/mod_bplot.R')

plot_perseus_volc <- function(df, group1, group2, curve.path,
                              y.max = NULL, x.min = NULL, x.max = NULL) {
  points <- df
  
  if(is.null(y.max)) {
    y.max <- max(-log10(points$pval))
  } 
  
  if(is.null(x.min)) {
    x.min <- min(points$Difference)
  }
  
  if(is.null(x.max)) {
    x.max <- max(points$Difference)
  }
  
  curve1 <- read_tsv(curve.path, show_col_types = F) |>
    filter(between(x, 0, x.max), y < y.max)
  
  curve2 <- read_tsv(curve.path, show_col_types = F) |>
    filter(between(x, x.min, 0), y < y.max)
  
  plabel <- points |>
    filter(type != 'ns')
  
  points  |>
    ggplot(aes(Difference, -log10(pval))) +
    geom_point(aes(color = type)) +
    geom_path(data = curve1, aes(x,y), alpha = .3) +
    geom_path(data = curve2, aes(x,y), alpha = .3) +
    geom_text_repel(data = plabel, aes(label = gene), color = 'purple') +
    theme_bw() +
    labs(title = str_c(group1, ' vs ', group2)) +
    scale_color_manual(values = c('lightgrey', 'red', 'blue'))
}

# SH2 domain protein -------
sh2.prot <- query_uniprot_keyword('KW-0727')

# PTD domain protein
ptd.prot <- c('prosite-PS51064', 'prosite-PS01179', 'prosite-PS51506') |>
  map(\(x)query_uniprot_keyword(x, qkey = 'xref')) |>
  bind_rows()

ptd.prot |> distinct(symbol)

# E3 protein ----------
e3.prot <- read_csv('mission/igg1_pY_omics/e3.ubiqitin.prot.csv')

# read xls -----------
perseus.volc <- 1:4 |>
  map(\(x)read_excel('mission/igg1_pY_omics/Perseus volcanos-IgG1-add MSMS.xlsx',
                     sheet = x))

## volcano with s0 curve ----------
logp_max <- max(c(perseus.volc[[1]]$...2, perseus.volc[[2]]$...2,
        perseus.volc[[3]]$...3, perseus.volc[[4]]$...3))

diff_range <- map(perseus.volc, pull, Difference) |>
  list_c() |>
  range()

perseus.volc[[1]] |>
  mutate(pval = 10^(-...2), Difference, gene = str_remove(`Gene names`, ';.+'),
         type = case_when(Significant == '+' & Difference > 0 ~ 'G396R.pY',
                          Significant == '+' & Difference < 0 ~ 'G396R.Y',
                          .default = 'ns') |>
           fct_relevel('ns'), .keep = 'none') |>
  plot_perseus_volc('G396R.pY', 'G396R.Y',
                    curve.path = 'mission/igg1_pY_omics/g2r.p.s0_2.txt',
                    y.max = logp_max, x.min = diff_range[1], x.max = diff_range[2])

g1 <- last_plot()

perseus.volc[[2]] |>
  mutate(pval = 10^(-...2), Difference, gene = str_remove(`Gene names`, ';.+'),
         type = case_when(Significant == '+' & Difference > 0 ~ 'WT.pY',
                          Significant == '+' & Difference < 0 ~ 'WT.Y',
                          .default = 'ns') |>
           fct_relevel('ns'), .keep = 'none') |>
  plot_perseus_volc('WT.pY', 'WT.Y',
                    curve.path = 'mission/igg1_pY_omics/wt.p.s02.curve.txt',
                    y.max = logp_max, x.min = diff_range[1], x.max = diff_range[2])

g2 <- last_plot()

perseus.volc[[3]] |>
  mutate(pval = 10^(-...3), Difference, gene = str_remove(`Gene names`, ';.+'),
         type = case_when(Significant == '+' & Difference > 0 ~ 'G396R.Y',
                          Significant == '+' & Difference < 0 ~ 'WT.Y',
                          .default = 'ns') |>
           fct_relevel('ns'), .keep = 'none') |>
  arrange(type) |>
  plot_perseus_volc('G396R.Y', 'WT.Y',
                    curve.path = 'mission/igg1_pY_omics/y.g2r.s02.curve.txt',
                    y.max = logp_max, x.min = diff_range[1], x.max = diff_range[2])

g3 <- last_plot()

perseus.volc[[4]] |>
  mutate(pval = 10^(-...3), Difference, gene = str_remove(`Gene names`, ';.+'),
         type = case_when(Significant == '+' & Difference > 0 ~ 'G396R.pY',
                          Significant == '+' & Difference < 0 ~ 'WT.pY',
                          .default = 'ns') |>
           fct_relevel('ns'), .keep = 'none') |>
  arrange(type) |>
  plot_perseus_volc('G396R.pY', 'WT.pY',
                    curve.path = 'mission/igg1_pY_omics/py.g2r.s02.curve.txt',
                    y.max = logp_max, x.min = diff_range[1], x.max = diff_range[2])

g4 <- last_plot()

g3 + g4 + g1 + g2

## volcano highlighting pY-binding protein ---------
perseus.volc[[1]] |>
  mutate(pval = 10^(-...2), Difference, gene = str_remove(`Gene names`, ';.+'),
         type = case_when(gene %in% sh2.prot$symbol ~ 'SH2 domain',
                          gene %in% ptd.prot$symbol ~ 'PTB domain',
                          .default = 'ns') |>
           fct_relevel('ns'), .keep = 'none') |>
  arrange(type) |>
  plot_perseus_volc('G396R.pY', 'G396R.Y',
                    curve.path = 'mission/igg1_pY_omics/g2r.p.s0_2.txt') +
  scale_color_manual(values = c('lightgrey','orange','green4'),
                     labels = c('other', 'PTB domain', 'SH2 domain'))

perseus.volc[[2]] |>
  mutate(pval = 10^(-...2), Difference, gene = str_remove(`Gene names`, ';.+'),
         type = case_when(gene %in% sh2.prot$symbol ~ 'SH2 domain',
                          gene %in% ptd.prot$symbol ~ 'PTB domain',
                          .default = 'ns') |>
           fct_relevel('ns'), .keep = 'none') |>
  arrange(type) |>
  plot_perseus_volc('WT.pY', 'WT.Y',
                    curve.path = 'mission/igg1_pY_omics/wt.p.s02.curve.txt') +
  scale_color_manual(values = c('lightgrey','orange','green4'),
                     labels = c('other', 'PTB domain', 'SH2 domain'))

perseus.volc[[3]] |>
  mutate(pval = 10^(-...3), Difference, gene = str_remove(`Gene names`, ';.+'),
         type = case_when(gene %in% sh2.prot$symbol ~ 'SH2 domain',
                          gene %in% ptd.prot$symbol ~ 'PTB domain',
                          .default = 'ns') |>
           fct_relevel('ns'), .keep = 'none') |>
  arrange(type) |>
  plot_perseus_volc('G396R.Y', 'WT.Y',
                    curve.path = 'mission/igg1_pY_omics/y.g2r.s02.curve.txt') +
  scale_color_manual(values = c('lightgrey','orange','green4'),
                     labels = c('other', 'PTB domain', 'SH2 domain'))

perseus.volc[[4]] |>
  mutate(pval = 10^(-...3), Difference, gene = str_remove(`Gene names`, ';.+'),
         type = case_when(gene %in% sh2.prot$symbol ~ 'SH2 domain',
                          gene %in% ptd.prot$symbol ~ 'PTB domain',
                          .default = 'ns') |>
           fct_relevel('ns'), .keep = 'none') |>
  arrange(type) |>
  plot_perseus_volc('G396R.pY', 'WT.pY',
                    curve.path = 'mission/igg1_pY_omics/py.g2r.s02.curve.txt') +
  scale_color_manual(values = c('lightgrey','orange','green4'),
                     labels = c('other', 'PTB domain', 'SH2 domain'))

g2r.p.sh2 <- perseus.volc[[4]] |>
  filter(`Gene names` %in% sh2.prot$symbol) |>
  mutate(log.pval = ...3, Difference, gene = str_remove(`Gene names`, ';.+'),
         .keep = 'none')

## volcano highlighting E3-associated protein ----------
perseus.volc[[1]] |>
  mutate(pval = 10^(-...2), Difference, gene = str_remove(`Gene names`, ';.+'),
         type = case_when(gene %in% e3.prot$gene ~ 'E3-associated',
                          .default = 'ns') |>
           fct_relevel('ns'), .keep = 'none') |>
  arrange(type) |>
  plot_perseus_volc('G396R.pY', 'G396R.Y',
                    curve.path = 'mission/igg1_pY_omics/g2r.p.s0_2.txt') +
  scale_color_manual(values = c('lightgrey','red'),
                     labels = c('other', 'E3-associated'))

g2r.p.sh2 |>
  mutate(gene = fct_reorder(gene, Difference),
         type = ifelse(Difference > 0, 'G396R-pY', 'WT-pY')) |>
  ggplot(aes(Difference, gene, size = log.pval, color = type)) +
  geom_point() +
  theme_bw() +
  geom_vline(xintercept = 0) +
  labs(title = 'SH2 containing protein: G396R-pY vs WT-pY',
       size = '-log10(p)')

perseus.volc[[4]] |>
  filter(`Gene names` == 'CRKL') |>
  select(contains('LFQ')) |>
  pivot_longer(everything()) |>
  mutate(group = str_extract(name, '(?<=_).+') |>
           fct_relevel('No_Probe')) |>
  ggplot(aes(group, value, color = group)) +
  stat_summary(geom = 'crossbar', fun = 'mean', fill = 'white') +
  geom_dotgraph +
  labs(y = 'LFQ intensity', title = 'CRKL') +
  theme_pubr()

g2r.p <- perseus.volc[[1]] |>
  mutate(pval = 10^(-...2), Difference, gene = str_remove(`Gene names`, ';.+'),
         type = case_when(Significant == '+' & Difference > 0 ~ 'g396r.pY',
                          Significant == '+' & Difference < 0 ~ 'g396r.Y',
                          .default = 'ns') |>
           fct_relevel('ns'), .keep = 'none')

gene.g2r.p <- g2r.p |>
  filter(type == 'g396r.pY') |>
  pull(gene)

wt.p <- perseus.volc[[2]] |>
  mutate(pval = 10^(-...2), Difference, gene = str_remove(`Gene names`, ';.+'),
         type = case_when(Significant == '+' & Difference > 0 ~ 'WT.pY',
                          Significant == '+' & Difference < 0 ~ 'WT.Y',
                          .default = 'ns') |>
           fct_relevel('ns'), .keep = 'none')

gene.wt.p <- wt.p |>
  filter(type == 'WT.pY') |>
  pull(gene)

# Venn of G396R-pY & WT-pY
list('WT-pY (vs WT-Y)' = gene.wt.p, 'G396R-pY (vs G396R-Y)' = gene.g2r.p) |>
  ggvenn(show_elements = T, label_sep = '\n', auto_scale = T,
         fill_color = c('blue','yellow'), text_size = 6)

py.g2r <- perseus.volc[[4]] |>
  mutate(pval = 10^(-...3), Difference, gene = str_remove(`Gene names`, ';.+'),
         type = case_when(Significant == '+' & Difference > 0 ~ 'G396R.pY',
                          Significant == '+' & Difference < 0 ~ 'WT.pY',
                          .default = 'ns') |>
           fct_relevel('ns'), .keep = 'none')

py.g2r.go <- py.g2r |>
  filter(!(gene %in% rna.splc) & type == 'G396R.pY') |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL',
           ont = 'ALL',readable = T, pvalueCutoff = 1, minGSSize = 3)

rna.splc <- py.g2r.go@result |>
  as_tibble() |>
  filter(Description == 'RNA splicing') |>
  pull(geneID) |>
  str_split_1('/')

py.g2r.go@result |>
  plot_enrichment(padj_thres = .1)

py.g2r |>
  filter(!(gene %in% rna.splc)) |>
  plot_perseus_volc(group1 = 'G396R-pY', group2 = 'WT-pY',
                    curve.path = 'mission/igg1_pY_omics/py.g2r.s02.curve.txt')


g2r.p.tidy <- g2r.pyvy |>
  select(`Gene names`, contains('G396R') & contains('LFQ'))

g2r.p.tidy |>
  pivot_longer(2:7) |>
  filter(`Gene names` == 'GRB2')

my.de <- g2r.p.tidy |>
  pivot_longer(2:7) |>
  mutate(group = ifelse(str_detect(name, 'pY'), 'pY', 'Y')) |>
  pivot_wider(names_from = group, values_from = value, values_fn = list,
              id_cols = `Gene names`) |>
  rowwise() |>
  mutate(t.test.p = t.test(pY, Y, var.equal = T)$p.value,
         avg_log2fc = log2(mean(pY) / mean(Y)),
         gene = `Gene names`, .keep = 'none')

g2r.pyvy |>
  mutate(pval = exp(-...2), Difference, gene = `Gene names`, .keep = 'none') |>
  left_join(my.de) |>
  ggplot(aes(pval, t.test.p)) + geom_point()

# SAM 10.1073/pnas.091062498
tusher.diff <- function(x, y, s0 = 0.1){
  nx <- length(x)
  ny <- length(y)
  a <- ((1/nx)+(1/ny))/(nx+ny-2)
  si <- sqrt(sum((x - mean(x))^2 + (y - mean(y))^2) * a)
  di <- (mean(x) - mean(y)) / (si + s0)
  message('number of sample x: ', nx)
  message('number of sample y: ', ny)
  message('gene-specific scatter si: ', signif(si))
  message('relative difference di: ', signif(di))
  di
}
tusher.diff(c(28.4,28.2,25.3), c(34.7,34.0,35.4), s0 = 1)
tusher.diff(c(28.4,28.2,25.3), c(28.4,28.2,25.3))
